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Abstract 

A key a priori information used in 4DVar is the knowledge of the system's evolution equations. In 
this paper we propose a method for taking full advantage of the knowledge of the system's dynamical 
instabilities in order to improve the quality of the analysis. We present an algorithm, four-dimensional 
variational assimilation in the unstable subspace (4DVar-AUS), that consists in confining in this 
subspace the increment of the control variable. The existence of an optimal subspace dimension 
for this confinement is hypothesized. Theoretical arguments in favor of the present approach are 
supported by numerical experiments in a simple perfect non-linear model scenario. It is found that 
the RMS analysis error is a function of the dimension A'' of the subspace where the analysis is confined 
and is minimum for N approximately equal to the dimension of the unstable and neutral manifold. 
For all assimilation windows, from 1 to 5 days, 4DVar-AUS performs better than standard 4DVar. 
In the presence of observational noise, the 4DVar solution, while being closer to the observations, 
if farther away from the truth. The implementation of 4DVar-AUS does not require the adjoint 
integration. 



1 Introduction 

Accuracy in the definition of the initial condition is an important factor for the performance of 
numerical weather and ocean prediction. The classical problem of estimating the state of a dynamical 
system from noisy and incomplete obse rvations is known in meteorology and oceanography as data 
assimilation (|Dalevl . 1 199]] : iKalnavl . [200^ . The goal of data assimilation in the initialization process is 
to provide the best possible estimate of the present state of the system using the available, partial and 
noisy, observations and the approximate equations governing the system's evolution. The estimate, 
referred to as the analysis, is obtained by optimally combining the information coming from a model 
forecast (background) and the observations (Talagrand, 1997). 

The non-linear stability properties of the system do not only determine the predictability horizon 
of the initial value problem but also profoundly infl uence the assim i lation process, affecting directly its 
quality and that of the subsequent forecast (see e.g. lCarrassi et al] . l2008al . and references therein). All 
assimilation methods, more or less implicitly, exert some control on the flow-dependent instabilities 
by means of the ob s ervat ional information. The Assimilation in the Unstable Subspace (AUS, 
iTrevisan and Uboldil |2004| ) explicitly estimates the flow-dependent instabilities and makes use of 
the unstable subspace as additional dy namical information. The 4- dimensional extension of AUS is 
the main scope of the present paper. In lTrevisan and Uboldil (|2004l ) and in other applications of the 
AUS assimilation, only a few unstable directions were tracked, whereas in the present study we make 
use of the entire unstable and neutral subspace, the subspace spanned by the Lyapunov vectors with 
positive and null exponents. 

Assimilation methods can be classified in two categories: sequenti al and variational, the most 
notable in the two c lasses being Kalman Filters and 4DVar respectively (|Ghil and Malanotte Rizzoli 
I1991I : iKalnavl . '2003', and references therein) . The Kalman Filter was originally developed for linear 
systems but a straightforward way o f extending the l inear results to th e nonlinear case is given by 
the Extended Kalman Filter (EKF) (|jazwinskil . 1l970l : iMiller et al.l . [l99i ') . 

Efficient minimization algorithms associated with adjoint techniques (|Talagrand and Courtier! . 

Il987l ) facilitate the implementation of 4DVar, an established and powerful assimilation method for 
meteorology and oceanography. In many realistic circumstances, reduced-rank approximations or a 



Monte Carlo approach, the latter referred to as Ensemble Kalman Filter (EnKF) (|Evensenl 1 1994 ). 
have been ad opted to circum vent the prohibitive cost of the full Extended Kalman Filter. The reader 
is referred to iKalnavl (|2003h and lTsuvuki and Mivoshil 12007) for a review o n the state of the art o f 
data assimilation in meteorology; see also Loren?l|2003 ) jKalnav et alj (|2007l ) and iGustafssonI l|2007h 
for a discussion on the relative merits of 4DVar and EnKF. 

The system's unstable subspace and its role in the assimilation process are central to our 
discussion. Hence, we briefly comment on how the flow dependent instabilities are dealt with 
by Kalman type fllters and 4DVar. In the Kalman filter, the propagation of the flow dependent 
instabilities is obtained by explicitly evolving the analysis error covariance from the previous analysis 
step. In Ensemble Kalman filters the subspace dimension of forecast error is at most equal to N-1, if N 
is the number of ensemble members: the rank deficienc y of the background error information is partly 
alleviated by covariance localization if N is too small (|Hamill et al.l . [2OOII I. A re lated aspect is the 
filter divergence that appears particularly critical in relation to sampling error (|Whitakcr and Hamillj, 
I2OO2I ). 

4DVar generates a model trajectory that best fits the observations available within a given 
assimilation window. Within the assimilation window, the fiow dependent instabilities are naturally 
described by the forward integration of the model and backward integration of the adjoint that model 
the error evolution. In addition, at the start of each ass imilation window, an a priori estimate of the 
background error covariance is needed l|Bannisteij , l2008t ) . 

For long assimilatio n windows, 4DVa r analysis errors are known to project on the unstable 
subspace of the system (jPires et al.l . Il996l ) . Errors in the stable directions that would be damped in 
the long range, for short assimilation windows are non- negligible i n the analysi s and affect the next 
assimilation cycle, causing short term enhanced error growth ( jSwanson et al ], iooo). It therefore 
seems appropriate to avoid introducing such type of error: this can be achieved by confining the 
increment of the 4DVar control variable in the unstable and neutral subspace of the system. In this 
way we avoid reintroducing observational error in the stable directions at each assimilation step: we 
anticipate that this is beneficial only if observations are not perfect. In this paper we present an 
algorithm (4DVar-AUS) that minimizes the 4DVar cost function under the above constraint. The 
dynamical information on the growth of errors in the unstable and neutral directions, the Lyapunov 
vectors with positive and null exponents, is explicitly estimated and, as explained in Section 12.2.21 
the adjoint integration is not needed. 

The idea of confining the analysis increment in the unstable subspace is not new . The sequential 
algorithm, referred to as AUS, has been introduced bv lTrevisan and Uboldil l|2004h . Its application 
to different models and observation configurations has shown that, even in the context of high- 
dimensional systems, an efficient error control and accurate state estimate can be obtained even by 
monitoringonly a reduced number of unstable directions (Uboldi and Trcvisan, 2006; Carrassi et al.„ 
[2OO7I . l2008bl ). The basic elements of the AUS scheme that differentiate it from other ensemble type 
Kalman filters are the explicit monitoring of the unstable directions of the system and the confinement 
of the analysis increment in the subspace that they span. A localization of the unstable structures, and 
consequently a localization of the error covariances (a feature common to other EnKF type methods) 
is necessary if the dimension of the subspace for the AUS assimilation is too small to describe the 
background error. The present extension of AUS to the four-dimensional case has the advantage of 
using the time distributed observations to track the instabilities that develop along the flow. 

One of the main goals of the present study is to address the following question: is there an optimal 
subspace dimension for the assimilation and is this related to the dimension of the unstable subspace 
of the system? In order to address this question, 4DVar-AUS is formulated in a perfect model setting 
and using a subspace of variable dimension. 

Theoretical arguments will be presented to indicate that the subspace dimension should at least 
be equal to the unstable manifold dimension. 

Results of the application of 4DVar-AUS to simple perfect non-li near sys tems obtained by varying 
the number of degrees of freedom in the Lorenz 40-variable model (|Lorenzl. [1996. 1 will be presented. 
The relation between the optimal subspace dimension and the number of positive exponents will be 
numerically investigated and it will be shown that the results conflrm the theoretical arguments. 

The paper is organized as follows: the formulation of 4DVar-AUS and theoretical argume nts on 
the subspace optimization are introduced in Section [2] results of the application to the lLorenzj l| 19961 ) 
model are presented in Section O while conclusions are drawn in Section [l] 



2 Formulation of 4DVar in the unstable subspace 
2.1 4DVar 

Strong constraint 4DVar seeks the (nonlinear) best estimate of the initial state xo that minimizes the 
misfit with observations in a given time interval (window) and possibly with a background state Xq. 
The standard cost function for strong constraint 4DVar, in discrete form, can be written as: 
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J(xo) = (xo - Xo)^B ^(xo-Xo) + 

+ f^(H,x.-yn^R"'(7i.x,-y°) 

i=0 

where y° are the observations available at discrete times ti = iAt, i = 0, n , within the assimilation 
window of length r = tn — to; B and R represent the background and observation error covariance 
matrices, TC the nonlinear observation operator, and the sequence of model states x^ is a solution of 
the nonlinear model equations: 

X, = Mo^^{xo), (2) 

The control variable for the minimization is the model state xo at the beginning of the assimilation 
window. 

Given the tangent linear equations describing the evolution of infinitesimal perturbations Sxi 
relative to an orbit of Eq. ((2)1 : 

5xj = Mo-,i5xo, (3) 
the gradient of J with respect to xq can be written as: 

^Vxo J = B~^(xo - Xo)+ 

(4) 

where represents the linearized observation operator, and the superscript T stands for 
transpose. 

For a given nonlinear trajec tory, the gradient can be estimated by use of the adjoint method 
(|Le Dimet and Talagraiidl . 119861 ). 

The solution of the minimization problem is obtained by forward integration of the model and 
backward integration of the adjoint, with an iterative descent algorithm. 



2.2 4DVar-AUS 

2.2.1 Unstable subspace and computation of Lyapunov vectors 

Lyapunov vectors, defined for nonlinear systems, are the time dependent physical structures 
associated with the Lyapunov exponents. There are basically two definitions of Lyapunov vectors 
that span the same invariant Oseledec subspaces. Th e first is that of an orthonormal set of vectors, 
the eigenvectors of the limit operator (jOseleded . ligGSl ): 

*ooW= lim [Mto->tMl^,]^(^ (5) 

to — * — oo 

where M* is the adjoint operator and the initial state, xo of the nonlinear trajectory is on the attractor 
(for fu rther reference in the meteorological literature see [Loronz (1984) and Lcgras and Vautarj 
\l99di)). 

The second is a set of non-orthogonal Lyapunov vecto rs, that are independ e nt of the norm and 
map into themselves with the tangent linear propagator l|Eckmann and RueUd . ITgSSI : iBrown et al.l . 
I1991I ) . These vectors hav e been shown to be the natural generalization of eigenvectors and Floquet 
vectors to aperiodic fiow (jTrevisan and Pancottil[l998l ). 

The following standard techn i que i s commonly used for calculating the orthonormal set of 
Lyapunov vectors l|Benettin et al.l . Il980h : a set of A'^ initially random tangent vectors are linearly 
evolved and orthonormalized every t time units. After a spin-up time these vectors span the A'^- 
dimensional most unstable subspace of the system. 

An efficient meth o d for recovering norm-independent non-orthogonal Lyapunov vectors is given by 
IWolfe and SamelsonI (|2007l ). Either one of the two above-mentioned methods can be used to identify 
the unstable, neutral, stable subspaces, the span of Lyapunov vectors with positive, null, negative 
exponents; the former, simpler techniqu e will be adopted in the prese nt application. 

In weather prediction, bred vectors l|Toth and Kalnavl [l993l . Il997l ') are usually computed instead 
of Lyapunov vectors. Bred vectors are the finite amplitude generalization of Lyapunov vectors 
and are computed as differences between twin nonlinear model integrations. The re-normalization 
amplitude and breeding time are the parameters being tuned to select the instability scale. With an 
infinitesimal re-normalization amplitude and periodic orthonormalizat ion the bred ve c tors a lgorithm 
would produce the same results as the Lyap unov vectors algorithm of lBenettin et all (|l980l '). 

In previous applications of AUS (see e.g. ICarrassi et al.l . l2008bl . and references therein) to more 
realistic atmospheric and oceanic models, the unstable vectors were computed with the breeding 
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technique. In those works, only a small number of bred vectors was used at each assimilation time. 
Because of the low dimensionality of the subspace spanned by those bred vectors a space localization 
was needed; furthermore, the breeding method was applied to the assimilation system instead of the 
freely evolving system in order to select those instabilities that survived the previous assimilation. 

In the present application we use the entire unstable and neutral subspace of the freely evolving 
system and we do not need any localization. 



2.2.2 The 4DVar-AUS algorithm 

The approach consists in determining the increment Sxo which minimizes the cost function in the 
reduced dimension subspace spanned by the A*' most unstable directions of the system corresponding 
to the leading A'^ Lyapunov exponents. 

After a transient time, the numerical technique described in Section 12.2.11 (|Benettin et al.l . Il980h 
leads to the identification of an orthonormal set of vectors spanning the A'-dimensional most unstable 
subspace. We apply this procedure to the solution of the assimilation cycle, starting initially with A'^ 
arbitrarily chosen tangent vectors. The Gram-Schmidt orthonormalization is applied at the end of 
each assimilation interval r. 

Let Eo be the matrix whose columns are the orthonormal tangent vectors spanning the A'- 
dimensional most unstable subspace of the system at to. The linear evolution within the assimilation 
window [tQ,T] is given by: 

Mo^^Eo = E,A,,(i = 0,...,n) (6) 

where 



exp 



Ai — diag 

,...,exp I X''"'>(t)dt 



to 
(iV) 



\^^Ht)dt,exp I \'^^'>{t)dt,..., 

(7) 



and X^-'\t) is the j"* local exponent. 

Consider the projection of the increment (Sxq in the subspace defined by Eq. In general, given 
a norm defined by the symmetric positive definite matrix Q, the projection can be written as 
Eo(Eo^QEo)^^Eo^Q(Sxo. For simplicity, in the following we adopt the Euclidean norm and we recall 
that the columns of Eo are orthonormal vectors. 

Thus, let the increment (5xo be confined in the subspace Eo and its projection (5xo be given by; 

S^o = EoE^'.Jxo (8) 
The evolution of the projected increment is governed by; 

= Mo_>EoE;[5xo = E, A.E^J'^xo, (i = 0, n) (9) 
Variations of the cost function JT]) due to variations of the control variable 5xo can be written as; 

SJ = (Vi^o J)^(5xo 
= (Vi^o J)'^(5xo 

where the tilde represents the projection into the subspace Eo, i.e. VxqJ = EqEq'VxoJ- 
Using Q and ([6]), the cost function gradient in the reduced subspace becomes; 



(10) 



-VkoJ" = Eo 



eJb-^(xo-x?,)+ 



(11) 



An assimilation cycle is constructed by initializing the next assimilation window with the state 
and its associated unstable subspace estimates at the end of the previous window. 
If k indicates the index of the assimilation cycle and recalling that t = tn — to: 
Xo(fc + 1) = XT(fc), where 

X,(fc) =Xo~.n(fc)(xo(fc)) (12) 

and Er(fc) = Eo(fc + 1)T, where 

Mo-.„(fc)Eo(fc) = E,(fc)A,(fc) (13) 

and T is the upper triangular orthonormalization matrix. 

In summary, the assimilation cycle is performed through the following steps; 
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1. A descent algorithm is used to find the cost function minimum by: 

(a) forward integration of the nonhnear model using ((2|) to compute Xi, starting from Xg at 
first iteration step; 

(b) forward integration of the perturbations Eo using ^ to compute and A^; 

(c) estimate of VxqJ from and of J from ([l]); 

2. The nonlinear model is integrated starting from the minimizing solution xo(fc) to produce the 
analysis, Xr(fc) H12|) . 

3. The perturbations Eo(fc) are evolved along the minimizing trajectory to produce ET-(fc) (|13|) : 

4. The columns of Er(fe) are orthonormalized and stored in Eo(fc + 1) to be used in the next 
assimilation cycle; 

5. Xo(fc + 1) is set equal to Xr(fc). 

Notice that no use of the adjoint integration is made. 

In the 4DVar-AUS assimilation, the analysis increment is confined in the A''- dimensional most 
unstable subspace of the previous analysis solution, with TV approximately equal to the number of 
positive and null Lyapunov exponents. 

Theoretical arguments given in Section [2]3l confirmed by numerical results presented in Section |3] 
will show that during a 4DVar-AUS assimilation cycle, errors in the stable directions are damped and 
errors in the analysis solution are confined withi n the unstable and neutral manifold of the system. 
This subspace is locally parallel to the attractor (|Eckmann and Ruellell 19851 ). so that one can find a 
state belonging to a nearby trajectory by moving along the tangent unstable directions. 

2.3 Full space and reduced order covariance matrix of the 
assimilation error 

Th e effect of the confi nement on the expected assimilation error covariance is now examined. 

iPires et al.l l| 19961 ) investigated the behavior of the observational term of the cost function in 
chaotic systems, making the tangent linear hypothesis and observing the whole state. They showed 
that, using the assumption that the observation error is uncorrelated in time and isotropic, with 
variance a^, the covariance matrix Co =< '7o^o"^ > of the assimilation error rjQ sd, t — to, <> being 
the expectation operator, can be written as: 



\i=0 / 

By confinement in the subspace defined by the N column vectors in Eo and using ((6)1 one easily 
obtains the following expression for the covariance of the assimilation error: 



To this point, no hypothesis has been made on the choice of Eo in (I15|l . If the number N of vectors of 
Eo is equal to the total number, /, of degrees of freedom of the system, (|15p represents the covariance 
matrix in the full space. 

Now let the N column vectors of Eq be the Lyapunov vectors corresponding to the N largest 
Lyapunov exponents. Assume that the Lyapunov vectors are orthogonal at to (or have been 
orthonormalized) and assume, for the sake of simplicity, that they remain orthogonal within the 
time window, then: 




(14) 



n 




(15) 



Co = CToEoDoEq' 




(16) 




where Ap' = expj^^ \^'\t)dt. 



At time t = t„, the analysis error covariance Ct =< VtVt'^ > under the tangent linear assumption. 



is: 



C, = 



cr^E^D.E^ 




(17) 




i=0 



5 



where Er = E(t„). 

In the expressions for Co and Cr, the role of the amphfying and decaying modes is interchanged. 
The generic term {Do)j.j of the diagonal matrix Do, representing the analysis error covariance 
associated with the j*'' Lyapunov direction (the j*'' column of E{to), ej{to)) at the beginning of 
the assimilation window, t = to is: 



+ ... + [exp{\^'^nAt)f} 

I _ g2(ii+l)A(3) At 



(18) 

2a(J> At 



where A'-*^ is the j*'' local Lyapunov exponent assumed, to simplify notation, to be constant within 
the assimilation window and At is the time interval between observations. A similar expression holds 
at the end of the assimilation window, t = t: 



{D.)j,j = + [exp{-\''''> At)f + [exp{-^\'^'hAt)f + 

+ ... + [exp{~\^'^nAt)f}~^ 

^ 2 1 - e 

° ]^ „ g-2(n+l)A(3) At 



(19) 



where, now, large n r efers to earlier tim es. 

In agreement with lPires et al.l (|l996l ). the influence of observational error on the stable (unstable) 
directions is damped as time increases (decreases) within the assimilation interval. At t = t 
(t = to), the largest error is along the most unstable (stable) directions; the older (more recent) 
observations, corresponding to increasing n, give smaller and smaller contributions. For sufficiently 
long assimilation windows the error along the stable (unstable) directions is damped. 



2.4 On the optimal subspace dimension 

Now, we focus on the effect on the analysis error covariance of the assimilation in the reduced, unstable 
and neutral subspace. In 4DVar-AUS, the analysis increment is confined in the subspace of the most 
unstable directions: let N be equal to A^ + A", where A^ and A" are the number of positive and 
null Lyapunov exponents. The influence of errors in stable directions, ejv+_|^jYO+i, e^+^jY0+2i --^i on 
the analysis error is eliminated by the confinement. Because we do not make corrections in the stable 
directions we avoid introducing, at every assimilation step, errors in the assimilation solution that 
project on the stable subspace. In this way, errors in the stable directions can naturally decay along 
the cycle. In standard 4DVar, instead, errors in the observations produce errors in the assimilation 
solution that project on the stable directions at each assimilation step. 

Results presented in the next section confirm these arguments, showing that the confinement in 
the unstable subspace is indeed beneficial to the performance of the assimilation. At final time of 
sufficiently long assimilation windows analysis errors in the stable directions would be damped also 
in standard 4DVar; however, this cannot be achieved in practice because nonlinearities set a limit to 
th e window extension th at can be used. 

ISwanson ^TZI dioooh investigated the effect of various observational errors on the 4DVar analysis 
and pointed out that errors in the stable directions can cause short term enhanced growth with 
adverse effect on the forecast. Their results are in agreement with the present findings. 

The effect of the confinement on the efficiency of the numerical algorithm is discussed in Appendix 



3 Application to the Lorenz (1996) model 

3.1 Model and experimental setting 

The results presented in th is section are based on the low-order chaotic model of iLoren^ (|l99(J l and 
iLorenz and Emanuell (|l998t ). The model is a simple analogue of mid-latitude atmospheric dynamics 
and its variables represent the values of a meteorological quantity at I equally spaced geographic sites 
on a periodic latitudinal domain. 
The governing equations are: 

^Xj = {xj+i - Xj-2)xj-i -Xj+F (20) 
with j = 1, I. 

Following iLorenz and Emanuell l|l998l) we set the external forcing F = 8, a value giving rise 
to chaotic behavior. In this paper we consider three model configurations with different numbers 
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1=40, N*=13 1=60, N*=19 1=80, N*=26 




Figure 1 : Time average RMS analysis error at t = r as a function of the subspace dimension N for three 
model configurations: 7=40, 60, 80. Different curves in the same panel refer to different assimilation 
windows from 1 to 5 days. The observation error standard deviation is cTq ~ 0.2. 



of degrees of freedom I. For I = 40, 60, 80 the three systems have 13, 19, 26 positive Lyapunov 
exponents, respectively. The doubling time associated to the leading Lyapunov exponent A^^-* is, in 
all three systems, approximately equal to 2 days if the system time unit corresponds to 5 days. 

Observing system simulation experiments are performed in a perfect model environment: a 
trajectory on the attractor of the system is assumed to represent the truth. Observations are created 
by adding to the true state Gaussian distribute d random errors wi th variance cTq. 

The observational network is the same as in iFertig et all |20o3). An observation is placed in one 
out of four grid points at each observation time. The frequency of observation is 1.5 hours, and the 
observed grid points are rotated so that, in a six hours interval, all grid points are observed once. 

An analysis cycle is set up with contiguous assimilation windows so that the initial time to in 
one window corresponds to the last time t„ = r of the previous one, as described in Section 12.2.21 
Experiments are performed using assimilation windows r — 1,...,5 corresponding to 1,...,5 days. 
In each ex periment th e anal ysis cycle consists of 5000 consecutive windows. A Conjugate Gradient 
algorithm ("P ress et all 1 19921 . Chap. 10) is used for the minimization of the cost function at each step 
of the algorithm of both 4DVar and 4DVar-AUS. 

The time mean, over an assimilation cycle of 5000 windows, of the RMS analysis error obtained 
with standard 4DVar is compared with that obtained by means of the 4DVar-AUS algorithm. The 
latter is applied using a variable number A'^ of directions in order to find the optimal subspace 
dimension for the confinement. 

3.2 Results 

Following the theoretical approach of IPires et al.l (|l996h . a first set of experiments is performed 
without background term in the cost function. The 4DVar-AUS algorithm described in Section [2.2.21 
requires the definition of the background (i. e. the estimate produced by the previous assimilation 
cycle) and of the subspace Eq representing instabilities at the end of the trajectory in the previous 
window. Therefore, even in the absence of an explicit background term in the cost function, the 
solution relative to one window is dependent on the solution of the previous one. For the first set of 
experiments an assimilation cycle is set up by successive minimizations of the observational part of 
the cost function only. We recall that the initial guess of the control variable at t = to is equal to 
the analysis at t = r of the previous assimilation and, in the 4DVar-AUS algorithm, the N column 
vectors of Eq are obtained by orthonormalizing the vectors E^. 

The results of 4DVar-AUS are compared with standard 4DVar. The RMS analysis error is 
computed at the end of each assimilation window and averaged over 5000 successive assimilation 
windows. Experiments are performed with the three systems with / = 40, 60 and 80 degrees of 
freedom; the number of null Lyapunov exponents of these systems is A" = 1 and the number of 
positive exponents is A^ = 13, 19, 26, respectively. 

3.2.1 Error dependence on the subspace dimension 

Figure [1] shows the mean RMS error as a function of the dimension A of the subspace Eq. Each panel 
refers to a different model configuration, / = 40, 60, 80. When N = I the error is that of standard 
4DVar (one can either set A = 7 in the 4DVar-AUS scheme or use the standard 4DVar algorithm: 
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Figure 2: Time average RMS analysis error at < = r as a function of the length of the assimilation 
window for three model configurations: 1=40, 60, 80. Different curves in the same panel refer to a 
different subspace dimension N of 4DVar-AUS and to standard 4DVar. CTq = 0.2. 



the results are the same within numerical accuracy). The value of ao is set to 0.2 (the 'climatological' 
standard deviation for the system is about 5.1). 

When A*' is smaller than the number of positive Lyapunov exponents, A*'^, the 4DVar-AUS 
algorithm does not converge or gives very poor results. When is increased above this threshold, the 
error abruptly decreases and then gradually increases again up to the value obtained with standard 
4DVar {N = I). Recalling that the number of positive global Lyapunov exponents for / = 40,60 
and 80 is 13, 19 and 26 respectively, the error minimum is obtained in all three model configurations 
for a value of A'' approximately equal to A'^ + Af". Because the value of local Lyapunov exponents 
fluctuates around the respective global value, even moderately decaying directions can be locally 
expanding and a number TV a few units larger than A'^+ + A''^ is needed. 

The most important result is that the minimum value of the error is obtained for an optimal 
subspace dimension which is very close to the number of positive and null Lyapunov exponents of 
the three (40, 60 and 80- variable) systems. 

Notice the internal consistency of the results of Fig. [T] the value of the average RMS error is 
virtually the same in the three model configurations. In fact, dynamically, the three models are 
equivalent and have the same instabilities, but a difi'erent number A'^"'' of unstable directions are 
present in proportion to the extension of the spatial domain. Because the observational configuration 
is the same, with a number of observations proportional to the domain size, and using the value of 
A'' appropriate for each system (given the respective value of A"^), we obtain the same accuracy of 
the analysis solution. 

Figure |2] displays the experimental data as a function of r, to illustrate the improvement obtained 
with the 4DVar-AUS scheme for different assimilation windows. Because the stable directions have 
a negative impact on the quality of the assimilation, the error appears to decrease by successively 
discarding a larger number of (the most stable) directions. We argued that errors in the stable 
directions do not affect the analysis for long enough (relative to the decay rate) assimilation windows. 
In agreement with this conjecture. Fig. [2] shows that the improvement obtained by using the optimal 
value of A'^, and thus eliminating the stable components of the error, is largest for the shortest 
assimilation windows: the largest improvement, a 30% reduction of the error with respect to classical 
4DVar is obtained for the smallest r corresponding to one day, while the improvement becomes less 
significant for larger r, and is about 20% for the five days window. 

The experiments were repeated using larger values of ao = 1 and \/2 and similar results were 
obtained (not shown) except that, as expected, the error scales with ao. It was not possible to 
complete all the experiments with a window of five days and the larger value of cto because the 
algorithm failed due to increased nonlinearity; in addition, the error minimum is shifted to slightly 
larger values of N: these were, for instance equal to 17, 25 and 30 for the 40, 60 and 80-variable 
systems respectively when the largest observational error value {ao = %/2) was used. This can be 
explained by noting that, when the error of the initial guess is larger, also the estimate of the unstable 
directions is less accurate at t = to; in such case a slightly larger subspace is needed. 



3.2.2 Stable and unstable error components within the assimilation window 

The time dependence of the error within the assimilation window, shown in Fig. |3] confirms the 
theoretical results of Section r2.3l and provides further insight on the reasons why 4DVar-AUS performs 
better than 4DVar. In the figure we show the 4DVar and 4DVar-AUS total error and the error 
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Figure 3: Time average RMS error within 1, 3, 5 days assimilation windows as a function of = i — r, 
with (Jo = .2, lO^'^ for the model configuration / — 40. Left panel: 4DVar. Right panel: 4DVar-AUS 
with N — 15. Solid lines refer to total assimilation error, dashed lines refer to the error component in 
the stable subspace eig, 649. 



projection on the stable directions e^^+^j^o^i, ...,ei, averaged over 5000 consecutive assimilation 
windows. The results shown were obtained with two values of CTo = 10~^ and .2, small enough that 
the analysis error scales with the observational error as predicted by the tangent linear theory of 
Section [Ol 

Figure [3] shows that, according to the theory, the 4DVar error is relatively larger in the stable 
subspace at initial time and in the unstable and neutral subspace at final time. Instead, in 4DVar- 
AUS errors are very small in the stable directions and project almost totally on the unstable and 
neutral subspace. 

The 4DVar-AUS assimilation error is smaller than the 4DVar error particularly for short 
assimilation windows. 

Because the search of the minimum of the 4DVar cost function is conducted in the entire phase 
space, its minimum cannot be larger than the minimum of the 4DVar-AUS cost function: this is 
confirmed by experimental evidence, the 4DVar cost function being typically a few percent smaller 
than the 4DVar-AUS cost function. We conclude that the 4DVar solution is closer to the observations 
but farther away from the real trajectory than the 4DVar-AUS solution. This is due to the fact that, 
in 4DVar, errors in the stable directions are "kept alive" by the observational error. In 4DVar-AUS, 
errors in the stable directions being never corrected, are naturally damped along the assimilation 
cycle: as a consequence, on average errors project only on the long term Lyapunov vectors contained 
in the matrix Eq. 

Results obtained by setting ao — show that the analysis error tends to zero in a time span 
that is shorter for standard 4DVar than for 4DVar-AUS; thus with perfect observations the full space 
4DVar assimilation performs better, strengthening our conclusions. It is worth mentioning that the 
evolution of Eq is a key factor for the performance of the assimilation: in fact, experiments performed 
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by using any number N < I of random directions show that the error is always larger than the error 
of the full space 4DVar assimilation (A'^ — I). 

3.2.3 Inclusion of the background term 

The 4DVar-AUS algorithm, in the absence of an explicit background term, amounts to assuming that 
the background error covariance matrix B is infinite in the unstable and neutral space, and in the 
stable subspace. The full 4DVar-AUS algorithm, still in the absence of an explicit background term, 
amounts to assuming that the matrix B is globally infinite. For completeness, a set of experiments 
was conducted with the inclusion of an explicit finite background term in the cost function. The 
static background error covariance matrix B was optimized for each assimilation window by the 
following iterative procedure. Starting from an initial guess, B is updated at each iteration step with 
the covariance of the difference between the forecast and true state, estimated from an assimilation 
cycle of 1000 consecutive windows. The process is repeated until convergence is obtained; in practice 
the iteration stops when the analysis error, averaged over the 1000 windows cycle, converges to an 
approximately constant value. To reduce the burden of computations B is optimized for each window 
only for standard 4DVar (A'^ = /); the same matrix is used for the 4DVar-AUS experiments with the 
same window: in this way 4DVar-AUS is penalized, since its results could only improve if we used a 
matrix B specifically optimized for each given subspace dimension N < I. 

Results are shown in Fig. |3J to be compared with Fig. [T] (for I = 40) and with Fig. |3] (for 
(To = 0.2). Everything else being equal, the introduction of the background term leads to an overall 
improvement of the performance in all experiments (compare the middle panel of Fig. |4] with the 
the lower left panel of Fig. [3] as concerns full 4DVar, and the right panel of Fig. |3]with the lower 
right panel of Fig. [3] as concerns 4DVar-AUS). This shows that useful information is contained in 
the matrix B. The most accurate analysis is still obtained for a value of A'' that is just above the 
number of positive and null Lyapunov exponents. It is seen that, in agreement with the theory, 
the largest improvements are obtained for shorter windows and for standard 4DVar. As expected, 
increasing the length of the time window decreases the influence of the background term. For 4DVar, 
the effect of the background term is particularly efficient in reducing the analysis error along the 
stable directions. Therefore a significant reduction of the error at the beginning of the time window 
is observed, but errors are nevertheless still present in the analysis in the stable directions. The 
conclusions to be drawn as to the matrix B as it has been defined here are first, as said, that it 
contains useful information in both the stable and unstable subspaces. Now, the error in the stable 
subspace is only partially decreased in the full 4DVar, while it is entirely eliminated in 4DVar-AUS. 
So, the description that the matrix B gives of the error in the stable subspace is more accurate than 
assuming that this error is infinite, but less accurate than assuming it is zero. 



4 Conclusions 

One of the main purposes of the present paper has been the development of four-dimensional 
variational assimilation in the unstable subspace. The results provide a proof-of-concept, at least 
in the case of the simple model used in this study, of the benefit in terms of assimilation performance 
of selecting the subspace where instabilities develop. The key result of this study is the existence of 
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an optimal subspace dimension for the assimilation that is directly related to the unstable and neutral 
subspace dimension. The selected subspace - the leading Lyapunov vectors subspace - contains the 
most rapidly growing perturbations. In the presence of observational error, the optimal number of 
directions is approximately equal to A^""" + A'"", where A'^"'" and A'^'^ are the number of positive and 
null Lyapunov exponents. The 4DVar solution (A'^ = 7), while being closer to the observations, is 
farther away from the truth. This result has been explained showing that, when we assimilate in 
the unstable and neutral subspace, errors in the stable directions are naturally damped. Because 
of observational error, assimilating in the whole space otherwise keeps the stable components of the 
error alive, deteriorating the overall assimilation performance. If the observational error is zero, the 
optimal dimension is the dimension / of the whole space. 

The present theoretical results can have implications for the application of advanced assimilation 
methods to high-dimensional models of the atmosphere and ocean. Here, we have shown that 4DVar 
could benefit from the dynamical information on the unstable directions - the "optimal" subspace 
where the analysis error is confined - propagated along the assimilation cycle. The possible application 
of the present findings to more realistic contexts is left for future investigations. Work is in progress 
to explore the existence of an optimal subspace dimension for EnKFs. 

The results presented in this paper have been obtained with a simple and economical numerical 
model and, strictly speaking, do not prove anything as to what would be obtained with a more realistic 
model of the atmospheric flow. At the same time, the considerations which have led to the design 
of the experiments described in the paper, and which are confirmed by the results that have been 
obtained, are very general, and do not fundamentally require anything else than the existence of both 
stable and unstable modes in the system?? under consideration. One possible source of difficulties 
could be the unavoidable presence of errors in the assimilating model. The experiments described 
here have been performed under the hypothesis of a perfect model. In the more realistic situation 
of an imperfect model , the corresponding e r rors w ill modify the unstable subspace, at least to some 
extent. The results of lTrevisan and Uboldil (|2004l l showed that the performance of AUS assimilation 
is not severely affected by the presence of model error in the model used above. Further work will be 
necessary in order to assess to which degree the presence of model errors in more realistic assimilation 
problems can affect the conclusions that have been obtained here. 

A Effect of the confinement on the efficiency of the 
numerical algorithm 

The minimization in the reduced subspace, chosen to be spanned by the most unstable directions is 
expected to converge more rapidly in view of the following argument. The Hessian, under the same 
simplifying hypothesis used to derive the analysis error covariance p4p can be written as: 



i=0 -' 

The condition number, ratio of the largest to the smallest eigenvalue of (|21|) . can be reduced if N 
is significantly smaller than the model space dimension. 
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